A New Method of Calculating the Spin- Wave Velocity c of Spin-1/2 Antiferromagnets 
With 0(N) Symmetry in a Monte Carlo Simulation 
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Motivated by the so-called cubical regime in magnon chiral perturbation theory, we propose 
a new method to calculate the low-energy constant, namely the spin-wave velocity c of spin-1/2 
antiferromagnets with O(N) symmetry in a Monte Carlo simulation. Specifically we suggest that c 
can be determined by c = L//3 when the squares of the spatial and temporal winding numbers are 
tuned to be the same in the Monte Carlo calculations. Here /3 and L are the inverse temperature 
and the box size used in the simulations when this condition is met. We verify the validity of 
this idea by simulating the quantum spin-1/2 XY model. The c obtained by using the squares of 
winding numbers is given by c = 1. 1348(5) Ja which is consistent with the known values of c in 
the literature. Unlike other conventional approaches, our new idea provides a direct method to 
measure c. Further, by simultaneously fitting our Monte Carlo data of susceptibilities xn an d spin 
susceptibilities \ to their theoretical predictions from magnon chiral perturbation theory, we find c 
is given by c = 1.1347(2) Ja which agrees with the one we obtain by the new method of using the 
squares of winding numbers. The low-energy constants magnetization density M and spin stiffenss 
p of quantum spin-1/2 XY model are determined as well and are given by M = 0.43561(l)/a 2 and 
p = 0.26974(5) J, respectively. Thanks to the prediction power of magnon chiral perturbation theory 
which puts a very restricted constraint among the low-energy constants for the model considered 
here, the accuracy of M we present in this study is much precise than previous Monte Carlo result. 
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INTRODUCTION 

During the last twenty years, models with 0{N) sym- 
metry which are relevant to antiferromagnets have drawn 
a lot of attention. In particular spin-1/2 Heisenbcrg-type 
models have been studied in great detail both analytically 
and numerically because it is believed that these models 
are the correct models to describe the undoped precursors 
of high T c cuprates (undoped antiferromagnets) . Beside 
their phcnomenological importance, these 0(N) models 
for antiferromagnets are interesting from theoretical per- 
spective as well. Further, because of the availability of 
efficient Monte Carlo algorithms and increasing comput- 
ing power, the physics of these models has been inves- 
tigated with unprecedented numerical accuracy fiTjlo|. 
For instance, using a loop algorithm as well as finite- 
volume and -temperature predictions from magnon chiral 
perturbation theory, the corresponding low-energy con- 
stants, namely the staggered magnetization density M s , 
the spin stiffness p s and the spin-wave velocity c of spin- 
1/2 Heisenberg model on the square lattices are deter- 
mined with high accuracy and are in agreement with ex- 
perimental results [ll|. Because the properties of these 
models with O(N) symmetry are well-studied, they are 
particular suitable for exploring any new idea. 

In analogy to chiral perturbation theory for the pions 
in QCD, a systematic low-energy effective field theory for 
the magnons in an antiferromagnet exists as well and is 
called magnon chiral perturbation theory [l2l4l4j . Low- 
energy effective field theories are based on symmetry con- 
straints of the underlying models and are universally ap- 
plicable. Results obtained by effective field theories are 



exact, order by order in a systematic low-energy expan- 
sion. Material specific properties enter the effective La- 
grangian in the form of a priori undetermined low-energy 
parameters, like the spin stiffness p (p s ) or the spin- wave 
velocity c. Once the numerical values of these low-energy 
constants are determined, either by Monte Carlo simula- 
tions or experimental data, the low-energy physics of the 
underlying models are completely determined. Since the 
low-energy physics of the underlying model only depends 
on the corresponding low-energy constants, it is impor- 
tant to determine these low-energy constants as precise 
as possible. From theoretical perspective, to determine 
the numerical values of these low-energy constants are 
important as well. For instance, by simulating spin-1/2 
Heisenberg model with an external staggered field on an 
exactly cubical space-time box (which requires a very pre- 
cise value of c), in addition to being able to determine 
the numerical values of A4 S and p s , such investigation 
also provides a good opportunity to exam the validity of 
predictions from the corresponding low-energy effective 
field theories fl5, 16|. 



For (sub-lattice) magnetization density and spin stiff- 
ness, one can directly measure the related observablcs 
and then use experimental finite lattice extrapolation for- 
mulae to obtain the bulk values of these 2 low-energy 
constants. On the other hand, the spin-wave velocity c 
is always determined in a less direct manner convention- 
ally. Motivated by the so-called cubical regime (defined 
later) in magnon chiral perturbation theory [13] , we pro- 
pose a new method to calculate the spin-wave velocity c 
of spin-1/2 antiferromagnets with O(N) symmetry in a 
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Monte Carlo simulation 1 . Specifically, we propose that 
c can be calculated by c = L//3 when the squares of 
spatial and temporal winding numbers are tuned to be 
the same. Here L and f3 are the spatial box size and 
the inverse temperature used in the simulations when 
above condition is met. Since this method allows one to 
measure c in a direct manner, the result is more accu- 
rate than other methods. Indeed as we will demonstrate 
later, for quantum spin-1/2 XY model, the numerical 
value of c = 1.1348(5) J a we obtain using the new idea is 
of high precision and is consistent with the known Monte 
Carlo results in the literature as well [lj|. Further, by 
simultaneously fitting our Monte Carlo data of suscepti- 
bilities xn and spin susceptibilities \ to their theoretical 
predictions from magnon chiral perturbation theory [l7| , 
we find c is given by c = 1.1347(2) Ja which is consistent 
with the one we obtain by the new method. These results 
confirm the validity and usefullness of our new method. 
Additionally other two low-energy constants, namely the 
magnetization density M. and spin stiffness p of spin- 
1/2 XY model are calculated with high accuracy and are 
much precise than previous Monte Carlo estimates [l9j ]. 

The remaining of this paper is organized as follows. Af- 
ter a brief introduction to our motivation of this study, we 
summarize the model and observables investigated here. 
Follows that we review the corresponding effective field 
theory predictions relevant to our study. Then we present 
our numerical results. In particular we demonstrate the 
validity of the method we used in our simulations to de- 
termine the low-energy constant c. Along the verifica- 
tion, the low-energy constants M. and p are also calcu- 
lated. Finally a section is devoted to the conclusion of 
our investigation. 



MICROSCOPIC MODELS AND 
CORRESPONDING OBSERVABLES 



The quantum XY model we consider in this study is 
defined by the Hamilton operator 



H 



(id) 



J 



~\~ Sjj s . 



(1) 



where Si and Sf are the first and second components 
of a spin- 1/2 operator at site i, and i, j denotes a pair 
of nearest neighbor sites on a square lattice. Further, J 
in cq. (Q]) is the antifcrromagnctic coupling. A physical 
quantity of central interest is the susceptibility which is 



1 This method was implicitly used in [IH for the study of con- 
straint effective potentials. Here we carry out quantitative in- 
vestigation to verified the validity of this method. 



given by 



1 



dt -Tr[M 1 (0)M 1 (i)exp(-/3 J ff)]. (2) 
Z 



Here /? is the inverse temperature, L is the spatial box 
size, Z = Trexp(— (3H) is the partition function and 
M 1 = £Y S x is the first component of magnetization. 
Another relevant quantity is the spin susceptibility which 
is given by 

1 f 13 1 

X = J dt -Tr[M 3 (0)M 3 (t)exp(-/iff)], (3) 

here M 3 = J2 X &x- Both xn and X can be measured very 
efficiently with the loop-cluster algorithm using improved 
estimators fill ]. In particular, in the multi-cluster version 
of the algorithm the susceptibility is given in terms of 
the cluster sizes \C\, i.e. \ii = ^77 (J2e l^l 2 )- Similarly, 
the spin susceptibility X = £ (wf) = £ (£ c W t (C) 2 ) 
is given in terms of the temporal winding number Wt = 
J2c Wt(C) which is the sum of winding numbers Wt(C) of 
the loop-clusters C around the Euclidean time direction. 
Finally, the spatial winding numbers are defined by Wi = 

£ c Wi(C) with »e {1,2}. 



LOW-ENERGY EFFECTIVE THEORY FOR 
MAGNONS 

Due to the spontaneous breaking of the global 0(2) 
symmetry, the low-energy physics of antiferromagnets 
with an 0(2) symmetry is governed by one massless Gold- 
stone boson. Detailed calculations of a variety of physical 
quantities for the spin-1/2 antiferromagnets with O(N) 
symmetry including the NNLO contributions have been 
carried out in [l7| . Here we only quote the results that 
are relevant to our study. The aspect ratio of a spatially 
quadratic space-time box with box size L is characterized 
by / = (/3C/L) 1 / 3 , with which one distinguishes cubical 
space-time volumes with (3c ~ L from cylindrical ones 
with /3c 3> L. The c appearing above is the low-energy 
constant spin-wave velocity. In the cubical regime, the 
volume- and temperature-dependence of the susceptibil- 
ity is given by 

M 2 L 2 p f c Q .,. 
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(4) 



where M. is the magnetization density and p is the spin 
stiffness. Further the spin susceptibility in the cubical 
regime takes the form 



(5) 
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In 



cq. ((4]), the functions /?,(/), which only depend 
on I, are shape coefficients of the space-time box de- 
fined in 17]. Finally, in the cylindrical regime, the 



temperature- and volume-dependence of the spin suscep- 
tibility x at very low temperature, namely when the con- 
dition L 2 p/{(3c 2 ) <C 1 is satisfied, is given by 



X 



2-^exp 



l c 2 /3 
2pL 2 



(6) 



THE DETERMINATION OF THE LOW-ENERGY 
CONSTANTS 

Conventionally the spin-wave velocity c is determi- 
nated indirectly in a Monte Carlo simulation. For ex- 
ample, magnon chiral perturbation theory predicts that 
the finite- volume dependence of the ground state internal 
energy density is given by [l7| 



e Q (L) =e + 0.7188725c/L 3 + 0(l/L 5 



(7) 



By fitting the related Monte Carlo data to above equa- 
tion, c can be obtained from the coefficient associated 
with the term of 1/L 3 in eq. (J7J|. Further, c can be cal- 
culated by the standard hydrodynamic relation \ = p/c 2 
as well. Notice for both methods mentioned above, ex- 
trapolations to infinite volume limit are necessary in or- 
der to obtain the numerical value of c. In other word, 
c is determined indirectly and the step of extraoplations 
will introduce systematic uncertainties into its numerical 
value. Because of this, here we propose a new method 
to calculate the numerical value of c for the spin- 1/2 an- 
tifcrromagnctic models with 0{N) symmetry directly in 
a Monte Carlo simulation. Our new idea is motivated 
by the so-called cubical regime in magon chiral pertur- 
bation theory [I?} . Specifically in magon chiral perturba- 
tion theory, an exactly cubical space-time box is obtained 
through the condition (3c = L, here again (3 and L are 
the inverse temperature and the spatial box size, respec- 
tively. Further, in a Monte Carlo calculation, if one sim- 
ulates the system in an exactly cubical space-time box, 
then the average of the square of spatial winding numbers 
(W 2 ) = \/2{{W'l) + (W£)) and the square of temporal 
winding number (W 2 ) should be the same. Once the 
condition (W 2 ) = (W 2 ) is met, c can be determined by 
c = L/f3. In practice, to employ this new method to cal- 
culate c, for a given box size L, one varied (3 until the con- 
dition (W 2 ) = (W 2 ) is reached. We would like to empha- 
size that the method we propose here to determine the 
low-energy constant c applies to any quantum spin-1/2 
antiferromagnetic system with a spontaneous symmetry 
breaking from a global O(N) symmetry to its 0(N — 1) 
subgroup. Notice since the c in the criterion of an exactly 
cubical space-time box, namely (3c = L is its bulk value, 
one would expect the c calculated by this new method 



suffers very mild finite lattice effects. Indeed, as we will 
demonstrate shortly, for the quantum XY model consid- 
ered in this study, th nuncrical value of c obtained by 
the new method is saturated to its bulk value even at 
L = 24a. 

To verify the validity of the new method we propose 
here to calculate c in Monte Carlo simulations, we have 
carried out several simulations using a continuum-time 
loop algorithm with L = 24a, 32a and L = 48a. Fur- 
ther, by tuning the inverse temperature (3 for each sim- 
ulations to reach an exactly cubical space-time box, the 
numerical values for c determined from these simulations 
with L = 24a, 32a, 48a, are given by c = 1.1349(8) Ja, 
c = 1.1347(10)Ja and c = 1.1347(7)Ja, respectively. 
Figure Q] demonstrates the results of such calculations. 
The 3 values for c obtained at different box sizes are 
consistent with each other and agree with earlier Monte 
Carlo result of c as well [ijj]. This provides a convinving 
evidence to support the validity of our new method of 
calculating c from the squares of spatial and temporal 
winding numbers. By a weighted average over these val- 
ues of c determined at different box sizes, the final re- 
sult of the numerical value for the low-energy constant 
c in this study obtained by the new method is given 
c = 1.1348(5) Ja. We have additionally carried out sim- 
ulations with L = 72a and (3 = 63.450099141/ J (which 
corresponds to c = 1.13475 Jo). The (W 2 ) and (W 2 ) 
obtained from these new runs are given by 17.1150(95) 
and 17.109(16), respectively. This result implies that the 
method of calculating c through the squares of spatial 
and temporal winding numbers indeed suffers very mild 
finite volume effects, at least for the model considered 
in this study. We notice that the p corresponding to 
these new runs is given by p = 0.26973(15)J which is 
statistically consistent with p = 0.26975(8) J calculated 
at L = 24a. 

Another way to exam whether the new idea of cal- 
culating c through the squares of winding numbers is 
quantitatively correct is to extract c by fitting Monte 
Carlo data of xn an d X i n the cubical regime to their 
predicted volume- and temperature-dependence, namely 
eqs. (j4} and j5]), respectively. However notice since c al- 
ways appears as the quantity c/p or p 2 /c in eqs. (j4]) and 
([5]), c and p are highly correlated and it would be a chal- 
lenge to extract c accurately from the fits. Fortunately 
we observe from our Monte Carlo data that the observ- 
able (W 2 ) which is exactly p(3 when L — > oo already 
reaches a constant for L > 24a. Hence in our simulta- 
neous fits, we also include an improved estimate of the 
(W 2 )//3 data obtained earlier when determining c using 
the new method and fit these data points to a constant. 
By simultaneously fitting cubical regime data of xn and 
X with L > 24a as well as the data of (W 2 )//3 to their 
predicted volume- and temperature-dependence formu- 
lae, we arrive at M = 0.43561(l)/a 2 , p = 0.26974(5) J 
and c = 1.1347(2)Ja with a x 2 /d.o.f. - 1. The results of 
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FIG. 1: The determination of c using the squares of spatial 
and temporal winding numbers. 



the fit are shown in figures and SJ The value of c cal- 
culated from the chiral fits agree nicely with the one de- 
tcmincd using the new method. Using large volume data 
points (L > 40a) for the fits leads to consistent results. 
This in turn proves the quantitative correctness of the 
new method of determining c using the squares of wind- 
ing numbers. Notice the values for M. and p we obtain 
are consistent with the known values from Monte Carlo 
simulations [l9j]. They are in good agreement with the 
related results from series expansion and spin- wave calcu- 
lations in the literature as well [2oT |. Further, the numer- 
ical values of these low-energy constants we obtain are 
much precise than those calculated in earlier Monte Carlo 




FIG. 2: Result of fitting the data points of xu obtained in 
the cubical regime to their magnon chiral perturbation theory 
prediction. Some data points are omitted for better visibility. 



study. Finally using the values for p and c determined in 
the cubical regime as well as eq. (|6|), we have compared 
the theoretical prediction and Monte Carlo data for \ m 
the cylindrical regime. The result of such comparsion is 
shown in figure [5] Considering the fact that there is no 
free parameter, the agreement demonstrated in figure [5] is 
reasonably good. All the results presented in this study 
also provides a strong support for the prediction power 
and quantitative correctness of magnon chiral perturba- 
tion theory in understanding the low-energy physics of 
the underlying model. 



CONCLUSIONS 

In this paper we have proposed a new method to cal- 
culate the low-energy constant, namely the spin- wave ve- 
locity c for general antiferromagnetic spin systems with 
a spontaneous symmetry breaking from a global O(N) 
symmetry to its 0(N — 1) subgroup using the squares 
of spatial and temporal winding numbers. We have 
demonstrated the validlity of this method by simulat- 
ing the quantum spin- 1/2 XY model. The numerical 
value of c we calculate with the new idea is given by 
c = 1. 1348 (5) J a which is consistent with the known 
Monte Carlo result in the literature. By fitting our 
Monte Carlo data of Xn and x to their volume- and 
temperature-dependence predictions from magnon chi- 
ral perturbation theory, we reach M. = 0.43561 (l)/a 2 , 
p = 0.26974(5)J and c = 1.1347(2)Jo. The value of 
c obtained from the fit is consistent with the one de- 
termined by the new method using winding numbers 
squared. This supports strongly the quantitative correct- 
ness of the new method we propose here to calculate the 
low-energy constant spin- wave velocity c in Monte Carlo 
simulations. The idea of using winding number squared 
is simple, but very powerful and requires moderate com- 
putational effort to obtain a very precise numerical value 
for c. Finally thanks to the robustness nature and predic- 
tion power of magnon chiral perturbation theory which 
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FIG. 3: Result of fitting the data points of (Wt) obtained in 
the cubical regime to their magnon chiral perturbation theory 
prediction. Some data points are omitted for better visibility. 



puts very restricted constraints on the low-energy con- 
stants and observables considered here, we are able to 
fit simultaneously our finite temperature data points to 
their predicted formulae from magnon chiral perturba- 
tion theory and obtain very accurate values for M, p 
and c. The agreement between theoretical prediction and 
Monte Carlo results of \ at very low temperature shown 
in figure O is remarkble as well considering the fact that 
there is no free parameter in obtaining figure [SJ 

ACKNOWLEDGEMENTS 

The simulations in this study were done based on the 
loop algorithms available in ALPS library [2l| ■ We would 
like to thank U. Gerber and W.-J. Wiesc for useful cor- 
respondence. Partial support from NCTS (North) is ac- 
knowledged. 



[I] B. B. Beard and U.-J. Wiese, Phys. Rev. Lett. 77 (1996) 
5130. 

[2] A. W. Sandvik, Phys. Rev. B 56, 11678 (1997). 

[3] B. B. Beard, R. J. Birgeneau, M. Greven, and U.-J. 

Wiese, Phys. Rev. Lett. 80 (1998) 1742. 
[4] A. W. Sandvik, Phys. Rev. Lett. 83, 3069 (1999). 
[5] Y. J. Kim and R. Birgeneau, Phys. Rev. B 62, 6378 

(2000). 

[6] L. Wang, K. S. D. Beach, and A. W. Sandvik, Phys. Rev. 

B 73, 014431 (2006). 
[7] F.-J. Jiang, F. Kampfer, M. Nyfeler, and W.-J. Wiese, 

Phys. Rev. B 78, 214406 (2008). 
[8] A. F. Albuquerque, M. Troyer, and J. Oitmaa, Phys. Rev. 

B 78, 132402 (2008). 
[9] S. Wenzel and W. Janke, Phys. Rev. B 79, 014410 (2009). 
[10] F.-J. Jiang, F. Kampfer, and M. Nyfeler, Phys. Rev. B 

80, 033104 (2009). 
[11] U.-J. Wiese and H.-P. Ying, Z. Phys. B 93, 147 (1994). 
[12] S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. 

Rev. B 39, 2344 (1989). 
[13] H. Neuberger and T. Ziman, Phys. Rev. B 39, 2608 
(1989). 

[14] P. Hasenfratz and F. Niedermayer, Phys. Lett. B268, 
231 (1991). 

[15] M. Gockeler and H. Leutwyler, Nucl. Phys. B350 (1991) 
228. 

[16] M. Gockeler and H. Leutwyler, Phys. Lett. B253 (1991) 
193. 

[17] P. Hasenfratz and F. Niedermayer, Z. Phys. B 92, 91 
(1993). 

[18] U. Gerber, C. P. Hofmann, F.-J. Jiang, M. Nyfeler, and 

U.-J. Wiese, JSTAT, P03021 (2009). 
[19] A. W. Sandvik and C. J. Hamer, Phys. Rev. B 60, 6588 

(1999). 

[20] C. J. Hamer, J. Oitmaa, and W.-H. Zheng, Phys. Rev. B 
43, 10789 (1991). 

[21] A. F. Albuquerque et. al, Journal of Magnetism and Mag- 
netic Material 310, 1187 (2007). 



p = 0.26974(5)J 



J 1 


\ 






♦ L=16 






• L = 20 






- L = 24 






• L = 28 








Ap = 0.26974(5)J : 






\c = 1.1347(2)Ja : 


-, T T 1 








10 20 30 40 50 60 70 



1000 1200 



FIG. 4: Result of fitting the data points of p obtained in the 
cubical regime to a constant. 



FIG. 5: Comparison between theoretical prediction and 
Monte Carlo data for \ m the cylindrical regime. 
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